/*!
 **************************************************************************
 * \file convdiff6.cpp
 * This code simulate Natural Convection in a unit cube with time-dependent
 * oscillating BC.
 *
 ***************************************************************************
 * \author    Luis M. de la Cruz [ Sat May 23 12:06:36 CDT 2009 ]
 ***************************************************************************
 * The equations to be solved are of the form:
 * \f[
 * \frac{\partial \phi}{\partial t} + \frac{\partial}{\partial x_j} 
 * \Big(u_j\phi\Big) =
 * \frac{\partial}{\partial x_j}\left(\Gamma\frac{\partial\phi}{\partial x_j}\right) + S,
 * \,\,\,\,\,\,\,\,\textrm{for  } j = 1,2,3,
 * \f]
 *
 * where \f$ \phi \f$ is a scalar variable \f$T, \rho, u_1, u_2\f$ or \f$u_3 \f$.
 * \f$ \Gamma \f$ is a diffusion coefficient and \f$ S \f$ is the source term.
 *
 * In this case we have an energy equation in terms of temperature \f$T\f$, 
 * coupled with the Navier-Stokes equations. 
 * The coupled equations are solved with the SIMPLEC strategy.
 * 
 * These equation are solved in a unit square
 * \f$ x, y \in [0,1] \times [0,1] \times [0,1]\f$. 
 * A fraction of two opposite walls have constant temperatures but different. 
 * The other walls are adiabatic. The no-slip condition is applied on all walls.
 * Next figure shows the domain and the oscillating boundary conditions imposed.
 *
 * \image html geom3Dtd.png "Oscillating boundary conditions" width=7cm
 * \image latex geom3Dtd.eps "Oscillating boundary conditions" width=7cm
 *
 *  To compile and run this example type the next commands:
 * \par
   \verbatim
   % make
   % ./convdiff6 \endverbatim
 * 
 ***************************************************************************
 */

#include <cmath>
#include <algorithm>

#include "Meshes/Uniform.hpp"
#include "Equations/ScalarEquation.hpp"
#include "Schemes/CDS_CoDi.hpp"
#include "Equations/Momentum_XCoDi.hpp"
#include "Schemes/CDS_XCoDi.hpp"
#include "Equations/Momentum_YCoDi.hpp"
#include "Schemes/CDS_YCoDi.hpp"
#include "Equations/Momentum_ZCoDi.hpp"
#include "Schemes/CDS_ZCoDi.hpp"
#include "Equations/PressureCorrection.hpp"
#include "Schemes/Simplec.hpp"
#include "Solvers/TDMA.hpp"
#include "Utils/inout.hpp" 
#include "Utils/num_utils.hpp"
#include "Utils/GNUplot.hpp"
#include "Utils/diff_operators.hpp"
using namespace Tuna;

// Scalar field in one dimension
typedef TunaArray<double,3>::huge ScalarField3D;

#define PI 3.14159265358979323846

//
// DT = 0.0001
//
//   int intervalo_1 = 50, intervalo_2 = 300;
//
// DT = 0.00001
//
//int intervalo_1 = 500, intervalo_2 = 3000;
//
// DT = 0.00001
//
//int intervalo_1 = 5000, intervalo_2 = 30000;

int intervalo_1, intervalo_2;
double freq;
int ciclo;

inline void boundaryConditions(ScalarField3D &T, int iteracion){

    int i,hot_up, hot_flat, hot_down;
    int it1 = intervalo_1;
    int it2 = it1 + intervalo_2;
    int it3 = it2 + intervalo_1;
    double thot, tcol;

    int cuarto_de_ciclo,octavo_de_ciclo,medio_ciclo;
    int op1,op2,op3;
    
    int num_vols_x = T.ubound(firstDim) + 1;
    int num_vols_z = T.ubound(thirdDim) + 1;
    Range I1(0,num_vols_x / 2 - 1), I2(num_vols_x / 2, num_vols_x-1);
    Range K1(0,num_vols_z / 2 - 1), K2(num_vols_z / 2, num_vols_z-1);

    hot_up = hot_flat = hot_down = 0;
    thot = tcol = 0.0;

    medio_ciclo = ciclo /2;
    cuarto_de_ciclo = ciclo /4;
    octavo_de_ciclo = ciclo/8;

    op1 = iteracion % ciclo;
    op2 = op1 % cuarto_de_ciclo;
    op3 = op2 % octavo_de_ciclo;

    /* up,flat,  or down*/
     if (0 < op3 and op3 <= it1){
          hot_up = 1;
          i = op3;
     }
     else if (it1 < op3 and op3 <= it2){
          hot_flat = 1;
          i=op3-it1;
     }
     else if (it2 < op3 and op3 <= it3 ){
          hot_down = 1;
          i=op3-it2;
     }
     
    /* thot or tcold */
    /* default: hot side is active */
    thot =  0.5 * ( hot_up * sin(freq*i) * sin(freq*i) +
		    hot_flat +
		    hot_down * sin(freq*(intervalo_1+i)) * 
		               sin(freq*(intervalo_1+i)) );

    tcol = 0;

     /* if cold side is active */
     if (octavo_de_ciclo < op2 and op2 <= cuarto_de_ciclo or op2 == 0) {
        tcol = -thot;
        thot = 0;
     }

     /* 1 2 3 or 4 sector */
     if (0 < op1 and op1 <= cuarto_de_ciclo){
	T(I2, T.lbound(secondDim),K1) = thot;
	T(I1, T.ubound(secondDim),K2) = tcol;
     }
     else if ( cuarto_de_ciclo < op1 and op1 <= medio_ciclo){
	T(I2, T.lbound(secondDim),K2) = thot;
        T(I1, T.ubound(secondDim),K1) = tcol;
     }
     else if ( medio_ciclo < op1 and op1 <= 3*cuarto_de_ciclo){
	T(I1, T.lbound(secondDim),K2) = thot;
	T(I2, T.ubound(secondDim),K1) = tcol;
     }
     else if ( 3*cuarto_de_ciclo < op1 and op1 <= ciclo or op1 == 0){
	T(I1, T.lbound(secondDim),K1) = thot;
        T(I2, T.ubound(secondDim),K2) = tcol;
     }
}

template<class T_energy, class T_xmom, class T_ymom, class T_zmom, class T_press>
double SIMPLEC(T_energy &, T_xmom &, T_ymom &, T_zmom &, T_press &, 
	       int, double);

double errorT, errorX, errorY, errorZ,errorP;

int main() 
{

    double lambda, Amplitud;
    double top_wall, bot_wall;
    double longitud_x, longitud_y,longitud_z;
    int num_nodos_x, num_nodos_y,num_nodos_z;
    double tolerancia, dt;    
    int max_iter, iteraciones_max, frecuencia;
    double alpha_p, alpha_u, alpha_v,alpha_w;
    double Prandtl, Rayleigh;

    int imon, jmon,kmon;

// Reading data from the input file
    std::ifstream input_cin ("input6");
    input_cin >> longitud_x
	      >> longitud_y
	      >> longitud_z
	      >> num_nodos_x
	      >> num_nodos_y
	      >> num_nodos_z
	      >> top_wall
	      >> bot_wall
	      >> dt
	      >> iteraciones_max
	      >> tolerancia
	      >> max_iter
	      >> Prandtl
	      >> Rayleigh
	      >> frecuencia
	      >> intervalo_1
	      >> intervalo_2;

    freq = PI / (2 * intervalo_1); 
    ciclo = 4*(4 * intervalo_1 + 2 * intervalo_2);
//
//  Linear System 
//
    SparseMatrix< Diagonal<double, 3> > A(num_nodos_x, num_nodos_y, num_nodos_z); 
    ScalarField3D             b(num_nodos_x, num_nodos_y, num_nodos_z);
//
//  Mesh definition
//
    StructuredMesh<Uniform<double, 3> > malla(longitud_x, num_nodos_x, 
					      longitud_y, num_nodos_y,
					      longitud_z, num_nodos_z);
    malla.print();
    double dx = malla.getDelta(X);
    double dy = malla.getDelta(Y);
    double dz = malla.getDelta(Z);
    int num_vols_x = malla.getExtentVolumes(X);
    int num_vols_y = malla.getExtentVolumes(Y);  
    int num_vols_z = malla.getExtentVolumes(Z);  
//
//  Scalar fields definitions (in volumes centers)
//
    ScalarField3D T(malla.getExtentVolumes());      // Temperature
    ScalarField3D p(malla.getExtentVolumes());      // pressure
    ScalarField3D pp(malla.getExtentVolumes());     // pressure correction
    ScalarField3D u(malla.getExtentVolumes());      // u-velocity
    ScalarField3D v(malla.getExtentVolumes());      // v-velocity
    ScalarField3D w(malla.getExtentVolumes());      // w-velocity
    ScalarField3D us(num_nodos_x, num_vols_y,num_vols_z); // u-velocity staggered
    ScalarField3D vs(num_vols_x, num_nodos_y,num_vols_z); // v-velocity staggered
    ScalarField3D ws(num_vols_x, num_vols_y,num_nodos_z); // w-velocity staggered
//
//  Scalar fields definitions for printing (in nodes)
//
    ScalarField3D Tn(malla.getExtentNodes());      // Temperature
    ScalarField3D pn(malla.getExtentNodes());      // pressure
    ScalarField3D un(malla.getExtentNodes());      // u-velocity
    ScalarField3D vn(malla.getExtentNodes());      // v-velocity
    ScalarField3D wn(malla.getExtentNodes());      // w-velocity
    ScalarField3D vorx(malla.getExtentNodes());   // vorticity
    ScalarField3D vory(malla.getExtentNodes());   // vorticity
    ScalarField3D vorz(malla.getExtentNodes());   // vorticity
//
// Energy equation definition
//
    ScalarEquation<CDS_CoDi<double,3> > energia(T, A, b, malla.getDeltas());
    energia.setGamma(1.0);
    energia.setDeltaTime(dt);
    energia.setNeumann(LEFT_WALL, 0.0);
    energia.setNeumann(RIGHT_WALL, 0.0);
    energia.setNeumann(FRONT_WALL, 0.0);
    energia.setNeumann(BACK_WALL, 0.0);
    energia.setDirichlet(TOP_WALL);
    energia.setDirichlet(BOTTOM_WALL);
    energia.setUvelocity(us);
    energia.setVvelocity(vs);
    energia.setWvelocity(ws);
    energia.print();
//
// Momentum-x equation definition
//
    Momentum_XCoDi<CDS_XCoDi<double, 3> > momento_x(us, A, b, malla.getDeltas());
    momento_x.setGamma(Prandtl);
    momento_x.setDeltaTime(dt);
    momento_x.setDirichlet(LEFT_WALL, 0);
    momento_x.setDirichlet(RIGHT_WALL, 0);
    momento_x.setDirichlet(TOP_WALL, 0);
    momento_x.setDirichlet(BOTTOM_WALL, 0);
    momento_x.setDirichlet(FRONT_WALL, 0);
    momento_x.setDirichlet(BACK_WALL, 0);
    momento_x.setUvelocity(us);
    momento_x.setVvelocity(vs);
    momento_x.setWvelocity(ws);
    momento_x.setPressure(p);
    momento_x.print();
//    
// Momentum-y equation definition
//
    Momentum_YCoDi<CDS_YCoDi<double, 3> > momento_y(vs, A, b, malla.getDeltas());
    momento_y.setGamma(Prandtl);
    momento_y.setDeltaTime(dt);
    momento_y.setDirichlet(LEFT_WALL, 0);
    momento_y.setDirichlet(RIGHT_WALL, 0);
    momento_y.setDirichlet(TOP_WALL, 0);
    momento_y.setDirichlet(BOTTOM_WALL, 0);
    momento_y.setDirichlet(FRONT_WALL, 0);
    momento_y.setDirichlet(BACK_WALL, 0);
    momento_y.setUvelocity(us);
    momento_y.setVvelocity(vs);
    momento_y.setWvelocity(ws);
    momento_y.setPressure(p);
    momento_y.setTemperature(T);
    momento_y.setRayleigh(Rayleigh);
    momento_y.print();
//    
// Momentum-z equation definition
//
    Momentum_ZCoDi<CDS_ZCoDi<double, 3> > momento_z(ws, A, b, malla.getDeltas());
    momento_z.setGamma(Prandtl);
    momento_z.setDeltaTime(dt);
    momento_z.setDirichlet(LEFT_WALL, 0);
    momento_z.setDirichlet(RIGHT_WALL, 0);
    momento_z.setDirichlet(TOP_WALL, 0);
    momento_z.setDirichlet(BOTTOM_WALL, 0);
    momento_z.setDirichlet(FRONT_WALL, 0);
    momento_z.setDirichlet(BACK_WALL, 0);
    momento_z.setUvelocity(us);
    momento_z.setVvelocity(vs);
    momento_z.setWvelocity(ws);
    momento_z.setPressure(p);
    momento_z.print();
//
// Pressure equation definition
//
    PressureCorrection<Simplec<double, 3> > presion(pp, A, b, malla.getDeltas());
    presion.applyBounds(1, num_nodos_x-1, 1, num_nodos_y-1, 1,num_nodos_z-1);
    presion.setUvelocity(us);
    presion.setVvelocity(vs);
    presion.setWvelocity(ws);
    presion.set_du(momento_x.get_du());
    presion.set_dv(momento_y.get_dv());
    presion.set_dw(momento_z.get_dw());
    presion.setPressure(p);
    presion.print();

    ofstream sorsumfile("./Data6/sorsum"), 
	nusselt_hot_file("./Data6/nuh"), nusselt_cold_file("./Data6/nuc"),
	umon_file("./Data6/umon"), vmon_file("./Data6/vmon"), 
        wmon_file("./Data6/wmon"), Tmon_file("./Data6/tmon"), 
        pmon_file("./Data6/pmon"),
	errorTfile("./Data6/errorT"), errorPfile("./Data6/errorP"),
	errorXfile("./Data6/errorX"), errorYfile("./Data6/errorY"),
        errorZfile("./Data6/errorZ");

    double sorsum, resorm, nuh, nuc;
    
    int iteracion;

    NumUtils::interpolateToNodes(Tn, T);
    NumUtils::interpolateToNodes(pn, p);
    NumUtils::interpolateToNodesUs(un, us);
    NumUtils::interpolateToNodesVs(vn, vs);
    NumUtils::interpolateToNodesWs(wn, ws);
    InOut::writeToFile_DX(pn, 0, "./Data6/pres.", dx, dy, dz);
    InOut::writeToFile_DX(T, 0, "./Data6/temp.", dx, dy, dz);
    InOut::writeToFile_DX(un, vn, wn, 0 , "./Data6/velc.", dx, dy, dz);
    
    for(iteracion = 1; iteracion <= iteraciones_max; ++iteracion) {
	boundaryConditions(T, iteracion);
	energia.updatePhi(T);
	
        sorsum = SIMPLEC(energia, momento_x, momento_y, momento_z, presion,
			 max_iter, tolerancia);

	if ( !(iteracion % frecuencia) ) {
	    NumUtils::interpolateToNodes(Tn, T);
	    NumUtils::interpolateToNodes(pn, p);
	    NumUtils::interpolateToNodesUs(un, us);
	    NumUtils::interpolateToNodesVs(vn, vs);
	    NumUtils::interpolateToNodesWs(wn, ws);

	    InOut::writeToFile_DX(Tn, iteracion, "./Data6/temp.", dx, dy, dz);
	    InOut::writeToFile_DX(pn, iteracion, "./Data6/pres.", dx, dy, dz);
	    InOut::writeToFile_DX(un, vn, wn, iteracion, "./Data6/velc.", dx, dy, dz);

	    DiffOperators::vorticity(vorx,vory,vorz, un, vn, wn, dx, dy, dz);
	    InOut::writeToFile_DX(vorx, vory, vorz, iteracion, "./Data6/vort.", dx,dy,dz);

	}

	if(!(iteracion % 1)) { 
	  
	  sorsumfile << iteracion << "\t" << sorsum << endl;
	  /***
	      nuh = NumUtils::calcNusselt(T, dx, dy, dz,0);
	      nusselt_hot_file << iteracion << "\t" << nuh << endl;
	  ***/

	  errorTfile << iteracion << "\t" << errorT << endl;
	  errorPfile << iteracion << "\t" << errorP << endl;
	  errorXfile << iteracion << "\t" << errorX << endl;
	  errorYfile << iteracion << "\t" << errorY << endl;
	  errorZfile << iteracion << "\t" << errorZ << endl;

	  umon_file << iteracion << "\t" << un(imon,jmon,kmon) << endl;
	  vmon_file << iteracion << "\t" << vn(imon,jmon,kmon) << endl;
	  wmon_file << iteracion << "\t" << wn(imon,jmon,kmon) << endl;
	  Tmon_file << iteracion << "\t" << Tn(imon,jmon,kmon) << endl;
	  pmon_file << iteracion << "\t" << pn(imon,jmon,kmon) << endl;
	  
	  cout << "=> Iter = " << iteracion 
	       << " s = " << sorsum 
	       << " nuc = " << nuc 
	       << " nuh = " << nuh << endl
	       << " e_T = " << errorT
	       << " e_P = " << errorP
	       << " e_X = " << errorX
	       << " e_Y = " << errorY << endl;
	}
    }

    std::cout << "\n\n NOTE: The data were written in the [./Data6] directory in DX format";

    return 0;
}

template<class T_energy, class T_xmom, class T_ymom, class T_zmom, class T_press>
double SIMPLEC(T_energy &energy, T_xmom &xmomentum, 
	       T_ymom &ymomentum, T_zmom &zmomentum, T_press &pressure, 
	       int max_iter, double tolerancia)
{
    double sorsum = 10.0, tolerance = 1e-2;
    int counter = 0;

    while ( (sorsum > tolerance) && (counter < 20) ) {
	energy.calcCoefficients();
	Solver::TDMA3D(energy, tolerancia, max_iter);
	errorT = energy.calcErrorL1();
	energy.update();

	xmomentum.calcCoefficients();
	Solver::TDMA3D(xmomentum, tolerancia, max_iter);
	errorX = xmomentum.calcErrorL1();
	xmomentum.update();

	ymomentum.calcCoefficients();
	Solver::TDMA3D(ymomentum, tolerancia, max_iter);
	errorY = ymomentum.calcErrorL1();
	ymomentum.update();

	zmomentum.calcCoefficients();
	Solver::TDMA3D(zmomentum, tolerancia, max_iter);
	errorZ = zmomentum.calcErrorL1();
	zmomentum.update();

	pressure.calcCoefficients();
	Solver::TDMA3D(pressure, tolerancia, max_iter);
	pressure.correction();
	sorsum = fabs( pressure.calcSorsum() );
	
	++counter;
    }

    return sorsum;
    
}

